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ABSTRACT 

Wc present the first non-LTE, co-moving frame molecular line calculations of a star- 
forming cluster simulated using smoothed particle hydrodynamics (SPH), from which 
we derive high-resolution synthetic observations. We have resampled a particle repre- 
sentation onto an adaptive mesh and self-consistently solved the equations of statistical 
equilibrium in the co-moving frame, using TORUS, a three-dimensional adaptive mesh 
refined (AMR) radiative transfer (RT) code. We verified the applicability of the code 
to the conditions of the SPH simulation by testing its output against other codes. We 
find that the level populations obtained for optically thick and thin scenarios closely 
match the ensemble average of the other codes. We have used the code to obtain 
non-LTE level populations of multiple molecular species throughout the cluster and 
have created three-dimensional velocity-resolved spatial maps of the emergent inten- 
sity. Line profiles of cores traced by N2H+ (1-0) are compared to probes of low density 
gas, -^^CO (1-0) and C^^O (1-0), surrounding the cores along the hue of sight. The 
relative differences of the line-centre velocities are shown to be small compared to the 
velocity dispersion, matching recent observations. We conclude that one cannot reject 
competitive accretion as a viable theory of star formation based on observed velocity 
profiles. 
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^ 1 INTRODUCTION 

K^ Large-scale hydrodynamic collapse calculations of turbulent 
^ clouds (e.g. Bate, Bonnell, & Bromm 2003 Krumholz, Klein, 



thin molecular tracers, to observe any motion of the ob- 



ject relative to the enveloping gas (Walsh, Myers, & Burton 



Cw |& McKee 2007 1 demonstrate the fragmentation of molecu 



20041 



lar clouds into clusters which in turn typically fragment into 
high-density cores and subsequently fragment into proto- 
stars. These simulations predict the statistical properties of 
the distribution of young stellar objects such as the initial 
mass function, the binary fraction and the star formation 
rate (e.g. |Zinnecker||1984| |Bonnell et al.||2003| |Bate et al 



Comparisons of the collapse calculations with observa- 
tions have generally relied on statistical tests of protostel- 
lar properties such as the IMF and binary fractions (|Bate| 



2009b but see also Kurosawa et al. 2004 1. Although un- 



doubtedly useful, such comparisons neglect the information 
encoded in the geometrical and dynamical properties of the 
star-forming gas and considering the wealth of observational 



2003[ |Bate||2009a| 7lt is expected that the fate of these ob- data available (e.g. |Narayanan et al.|2008| [Buckle et al.|2010 1 



jects is strongly affected by their evolution history, specifi- 
cally the competitive accretion of gas as they move through 



their environs. Star formation in isolation ( Adams, Lada, & 
Shu|1987l necessarily neglects the dynamic interactions be- 



tween these objects, which can cause them to be ejected at 
great velocities from their nascent gas-rich core, truncating 
their accretion (Klessen, Burkert, & Bate||T998 Bate, Bon- 



nell, fc Bromm|2002 I and it will be possible, using optically 



these tests should provide a further, robust appraisal of the 
models. In order to make comparisons between theory and 
observation, given the intractability of the inverse problem, 
it is necessary to transform the temperature and density 
characteristics of a simulated cloud into a synthetic intensity 
map comparable with those obtained by observation. That 
being so, this paper aims to improve upon the analysis of the 
relative motions of cores to their nascent envelopes in the 
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clust er models of|Bate et al.| ( |2003[) by|Ayliffeet al.| ( |2007 1 . 
Both [Aylifle et al.| ( |2007| ) and [Walsh et al.| ( |2004| | look for 
evidence of supersonic ballistic motions of cores relative to 
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their inchoate envelopes. Assuming that the velocities of the 
core and the envelope are in the line of sight it is possible 
to determine the difference between them. By performing 
a non-LTE radiative transfer calculation of the molecular 



level populations within the cluster using TORUS (see Harries 
2000[ [Harries etal. 120041 [Symington et al.|2005| for details), 



we seek to improve upon the comparison between simula- 



tion and observation carried out by Ayliffe et al. ( 2007 1 and 



reconcile the differences between observation and theory. 

In order to understand the chemical and physical con- 
ditions in these structures it is necessary to understand the 
role that molecules play in radiation transfer in these cool, 
dense regions. For over 40 years (e.g. |Dieter[Pl964[ |Wilson| 
et al.|1970 1 molecular line radiation has regularly been used 
to trace key physical parameters such as density, temper- 
ature and velocity. Today's (sub-) millimetre observatories 
enable the observation of some of the coolest, most chemi- 
cally rich regions of the galaxy. Furthermore, the Atacama 
Large Millimeter Array (ALMA) promises unprecedented 
resolution of these objects therefore it is vital that the in- 
terpretation of the data is sound, necessitating the use of 
non-symmetric models incorporating complex physics. This 



was first recognized by Bernes ( 1979 1 and today many line 
radiative transfer codes exist (e.g. Juvela|1997 Wiesemeyer 



|1997[|Rawlings fc Yates|2001| ). Until relatively recently, only 
spherically or axially symmetric core collapse sinmlations of 



line transfer had been performed (e.g. Tsamis ct al._ 2008( 
[Pavlyuchenkov et al. 12008] ) however, due to recent advances 
in computational power, it is now possible to compute three- 
dimensional line transfer for these complex density struc- 
tures in star-forming regions. Here we use the TORUS code, 
since the use of adaptive mesh refinement (AMR) enables us 
to treat the huge dynamic range in density and linear scale 
associated with the SPH cluster calculations. 

In section [2] we present our implementation of a non- 
LTE molecular line radiative transfer module for TORUS. 
In section [3] we demonstrate agreement with other codes 
in some of the benchmark problems set out in |van Zadel-| 
hoff et ah] ( [2002 1 and verify the accuracy of our raytracing 
routine for generating datacubes. In section [4] we present 
an efficient method for mapping irregular SPH data on to 
an AMR grid and in section [5] we present the results of its 
application to the hydrodynamic clustered star-formation 
simulation of Bate et al. (|2003[). Finally, we discuss the re- 



sults of our analysis of the cluster simulation and conclude 
on our findings and on the importance of radiative transfer 
(RT) analysis of star formation as a whole. 



2 METHOD 

In order to accurately determine the local radiation field 
due to the presence of molecules, our radiative transfer 
code adopts the accelerated Monte-Carlo (AMC) method 
described in |Hogerheijde fc van der Tak| ( |200()| ). The greatest 
benefit of AMC methods is the ability to solve the equations 
of radiative transfer in regions of exceptionally high opti- 
cal depth (r > 100). In particular, [Hogerheijde fc van der| 
|Tak[ s method uses rays, or long characteristics, to sample 
the external radiation field as opposed to tracking photon 
packets, c.f. |Lucy| ( |1999| ), which can become trapped lead- 
ing to poor convergence or sometimes failure to converge at 



all. Using this method, each grid cell is sampled individually 
meaning that no cell is undersampled. Furthermore, the sep- 
aration of local and external contributions to the radiation 
field in a cell facilitates convergence even in optically thick 
regions. This is analogous to the operator-splitting step in 



Accelerated Lambda Iteration (ALI) methods (see Rybicki 
fc Hummer||1991 for details) . 



2.1 Determining non-LTE level populations 

In the scheme, there are two stages; the first stage samples 
the intensity of radiation incident upon each cell system- 
atically using a small set of rays (we use 100) each pos- 
sessing a unique origin within the cell, direction (d) and 
a frequency (i/). The frequency is sampled from a uniform 
distribution with a width of 4.3 normalized turbulent line 
widths, vt-arhVo/c, Centred on the rest frequency for the tran- 
sition, vo. Beyond 2.15 linewidths, the line profile function, 
4>(v), drops to less than 1 per cent of its peak value. The 
distribution from which the frequency is picked is uniform 
to ensure good sampling in the optically-thin line-wings. In 
the first stage, the direction and frequency of the set of rays 
do not vary from iteration to iteration, nullifying the ran- 
dom fluctuations in coverage of the radiation field and conse- 
quently converges rapidly towards a self-consistent solution. 
However, because the radiation field is poorly sampled both 
spatially and in frequency, the solution will be strongly de- 
pendent upon the choice of initial random seed. 

The second stage of the scheme reduces these system- 
atic and random errors by doubling the number of rays used 
to sample the radiation field per cell per transition with each 
successive iteration. The noise in a solution from a Monte- 
Carlo simulation diminishes proportional to \/iV asymptot- 
ically, where A^ is the number of rays used to sample the 
field, so it is quite possible to have to sample a cell with in 
excess of 10® rays in order to reach an acceptable level of 
convergence. Moreover, by allowing the properties of each 
ray to vary from iteration to iteration any source of system- 
atic error will also be reduced. 



Juvela (19971 describes a number of schemes for reduc- 



ing the noise and accelerating the convergence of Monte- 
Carlo methods. We have used low-discrepancy sequences 
to optimally cover the direction/frequency parameter space 
and have implemented Ng acceleration (Ng 19741 with 



weights from Olson, Auer, fc Buchler ( 1986 \ to improve the 



convergence properties of our code without loss of generality. 
We have not yet implemented any kind of spatial weighting 
to reduce the variance between iterations for the code and 
as such it is possible that the contribution of a small bright 
grid cell to another cell at some distance will often be missed 
when the number of rays used per cell to sample the field 
is low. However, in the scenarios considered in this work, it 
has not been necessary to implement such a strategy. 



2.1.1 Determining the external radiation field 

Common to both stages is the calculation of the mean in- 
tensity, ,Jv, experienced by each cell in the grid. This is de- 
termined by averaging the contribution of samples of the 
specific intensity incident at a point in the cell over 47r stera- 
dians. The specific intensity, I^, upon the point of origin of 
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the ray then follows from integrating the radiative transfer 
equation, 

where 



J-u \ Ot/, 



(1) 



dTv — Uuds 



and 



S. = 



J^ 



(2) 



(3) 



in a piecewise manner along the ray direction towards the 
edge of the grid, the contribution of each subsequent cell 
being attenuated by the frequency specific optical depth, r^, 
of all intervening cells. The source function, Sv, is defined 
as the ratio of the emission coefficient, j^ (erg s~^ cm""^ 
Hz~^ sr~^), to the absorption coefficient, a^ (cm~^). We 
assume the blackbody radiation from the cosmic microwave 
background with characteristic temperature Tcmb = 2.73 K, 
B{TcMS, vq) as an external boundary condition, where B is 
the Planck function. 

All physical quantities within the cell are assumed to 
be constant except for the systematic velocity vector, v(x), 
which enters equations ( |2|3| l through the line function 4>{'^) 
in both the emission and absorption coefficients (in the ab- 
sence of dust) : 



0!„ — — — (mBiu — nuBui)4>v 
47r 



(4) 
(5) 



Aui, Bui and -B;„ are the Einstein coefficients of sponta- 
neous emission, stimulated emission and stimulated absorp- 
tion associated with the transition from level u to level I 
respectively; similarly Wi, and ni are the respective number 
densities of molecules populating those levels. Absorption 
due to dust can be parameterised by k,^(T) which may also 
vary as a function of position. Reprocessed dust emission at 
the same frequency is determined by the Planck function at 
the dust temperature, Tdust which is assumed to be equal to 
the gas temperature unless otherwise stated. 

The line profile function (?!>(i^), centred on frequency vq, 
is characterized by the line broadening parameter «turb: 



UturbfOV^ 



exp 



Av' 



VtVLTh = \ V£,+V 



(6) 



(7) 



recalling that the thermal broadening ^t ~ ^ where m is 
the molecular mass and where i^^t is used to parameterize 
the unresolved non-thermal component of the gas velocity. 
Av is the equivalent velocity required to Doppler-shift i^o 
to v. However the presence of a systematic velocity field 
introduces anisotropy by its inclusion in Aw, thus: 



Av = c- 



1^0 



vo 



+ V- d. 



(8) 



TORUS accepts an analytical function for evaluation of 
v(x). However, in the absence of one, or if the velocity 
is computationally expensive to calculate, it is possible to 
quadratically interpolate over values stored at the grid cell 



corners. The local source function can vary strongly with 
velocity gradient across a cell when the velocity gradient is 
large compared to the local magnitude of the turbulent ve- 
locity field. Variations in 4>{u) are determined by subdividing 
the integration into smaller steps. To ensure good resolution 
of the velocity structure and consequently the differential 
optical depth within a cell, we use the condition. 



Wapiit = max I 2, [5 



fturb 



(9) 



where Ve and Vs are the velocities at the cell boundaries 
where the ray intersects with the cell. 

Equationlllis integrated along each ray from cell bound- 
ary to cell boundary to the grid edge for each transition be- 
ing considered to build up a representation of the radiation 
field external to the cell (as well as quantifying the extent 
to which it has been attenuated by intervening material), 
which when added to contributions from within the cell ap- 
proximates the mean radiation field. 



2.1.2 Determining J and solving the equations of 
statistical equilihnum 

The mean radiation field, ,]„, for the cell is determined by 
averaging the contribution from each ray, weighted by its 
line profile function. Ji, is split into an external component 
whose value is fixed and an internal component that varies 
according to the relative ratios of the level populations which 
affect Sij so that 



fJh' — t'l/ \ <J II 



Ei^e- 



E, 



- + 



E.g.[l-e" 
E.<^. 



(10) 



In order to obtain the instantaneous relative fractional 
level populations, n^, for a cell it is necessary to solve a sys- 
tem of equations describing the statistical equilibrium that 
is to be attained, the so-called equations of detailed balance. 



ni 



'^Aik + ^i.BikJ^ + Cik 



ki^l 



k>l kjtl 



(11) 



where Cik are the total rate coefficients of collision from level 
I to k for the molecule and a particular collision partner (e.g. 
H2 or e~) at a specific temperature. 

The results of non-LTE calculations are known to 
be highly sensitive to the quality of these coefficients 
so good estimates are essential. TORUS uses molecular 
data {vo, ^ui, Bui, Bia,Ciu{T)) from the LAMDA database 
( Schoier et al.|[2005 l, a repository of molecular and atomic 
data for many common species. 



For each cell, the system of equations (10 111 are solved 



iteratively to obtain a self-consistent solution vector of rela- 
tive level populations. TORUS uses a user-defined parameter 
to specify the point at which successive iterations are de- 
termined to have converged. By default, the criterion is set 
so that the root-mean-square (RMS) error between the two 
most recent iterations over all but the two uppermost levels 
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is less than 10 ^^ , i.e. 



\ 



E 



< 10" 



(12) 



Typically, the uppermost levels are chosen so that their pop- 
ulation has no significant impact on the quality of the solu- 
tion of lower levels. Due to their low population, the upper- 
most levels are often subject to a high degree of noise and 
their convergence is not representative of the convergence of 
the other levels. This criterion ensures every important level 
is converged to at least 10"^" and many will be converged 
to a much higher degree. 



2.1.3 Convergence and acceleration 

Once each cell has an updated set of level populations, they 
are compared with those from the previous iteration. The 
grid is considered to be converged when the RMS error over 
all cells is less than a user-specified tolerance, typically 1 per 
cent, for all levels. 

Once the grid has converged using fixed sets of rays, the 
algorithm enters the second stage using increasingly large 
sets of rays to sample the mean radiation field. Each iter- 
ation may potentially take a long time to complete; conse- 
quently, at the end of each iteration the intermediate grid is 
saved along with diagnostic information that can be useful 
to determine if there are pockets of slower convergence in an 
otherwise converged grid. 

As with many problems that are solved numerically by 
repeated function evaluation it is often possible to obtain 
a more accurate solution more rapidly by employing vector 
sequence acceleration. TORUS uses a vector acceleration tech- 



nique developed by Ng ( 1974 1 to extrapolate an updated set 



of relative level populations from the previous four iterations 
of rii 



Olson et al. 



( 1986 ) comment that weighting each level 
with Wi = Jul ~ ensures regions with small J are adequately 
represented. A further advantage of this method is that be- 
cause the vector is taken from a space where ^rii = 1, 
the relative level populations will still be normalized. In a 
complex iterative scheme such as this one where function 
evaluation can be expensive, we find that these acceleration 
techniques, used appropriately, can not only reduce the time 
required to find an accurate solution but may, in cases of ex- 
tremely high optical depth, be used to determine a solution 
that might not be possible to obtain by fixed iteration at 
all. 



2.2 Datacube generation 

As well as determining the non-LTE level populations of 
a simulation, TORUS can create three-dimensional velocity- 
resolved spatial maps of the emergent intensity (hereinafter 
described as datacuhes). Datacubes of any size may be cre- 
ated, limited only by the amount of memory and the spatial 
resolution of the original AMR grid. Additionally, TORUS 
stores extra information such as the optical depth reached 
in a particular datacube element and the column density. 
TORUS generates datacubes in a similar way to that in which 
it finds a non-LTE solution for the level populations within 
a grid. The incident radiation at a point is determined by 



querying the source function along the path of the ray, ex- 
cept that when creating datacubes the point lies outside 
the grid, on an image plane centred on the hypothetical 
observer's position and that now the frequency of the ob- 
servation is chosen by the observer. Typically, the datacube 
is comprised of many velocity channels which the observer 
may scan through in order to create images of intensity at 



each corresponding frequency. Douglas et al. (20101 have 



extended this code to render images in solid angle (as op- 
posed to rectilinear coordinates) in the near field where the 
observer is close to or actually inside the object, that is, ob- 
serving the spiral structure of a galaxy in Hi from within 
the galaxy itself. 

The plane of projection is defined by the vector passing 
through the point of observation and a point inside the grid. 
Orthonormal basis vectors for the image are chosen such 
that one is perpendicular to z (in the case where the observer 
is looking along z, the bases are chosen as ±x and ±J/). The 
image is then subdivided into A''^ square pixels of equal area 
so that the intensity per pixel may be stored in the datacube 
structure. 

More often than not, the resolution of the image will 
be less than that of the grid being imaged, i.e. the width 
of a pixel will be greater than that of the smallest subcells 
contained within the grid. This necessarily means that some 
information about the structure will be lost. This is espe- 
cially likely to occur if the intensity in a pixel is assumed 
to be that calculated by the pencil-beam emanating from 
the centre of the pixel integrated over the area of the pixel. 
TORUS utilizes sub-pixel sampling to ameliorate this effect, 
by keeping a running average of the intensity sampled over 
different paraxial rays bound within the pixel until the stan- 
dard error is less than some user-specified amount (we use 1 
per cent) or until the number of rays exceeds a pre-specified 
maximum. 

The origin of the ray inside the pixel is chosen using a 
Sobol quasi-random number generator that is reset for each 
new pixel. The advantages over pseudo-randomly generated 
origins are two-fold; intuitively, the first ray always origi- 
nates at the pixel centre, and, the origins are created so 
that they avoid each other as best possible within the con- 
fines of the pixel thereby avoiding dumpiness (mathemati- 
cally, a property known as low discrepancy) . This provides a 
truer representation of the range of intensities in the region 
being sampled. Furthermore, because these low discrepancy 
sequences are capable of being extended in a pointwise fash- 
ion until some criterion is met, they are considerably more 
flexible than grid-based subsampling techniques that require 
an exponential increase of samples. Moreover, in geometries 
where the density can be evaluated analytically or is stored 
on the corners of each cell in addition to the cell centres, 
TORUS uses density subsampling. 

In the iterative level population solver, a ray is subdi- 
vided within a cell to ensure good resolution of the velocity 
field. Analogously, using this technique it is possible to use a 
more accurate value for the density at a point along the ray 
than that of the cell centre itself producing both smoother 
and more accurate images at lower grid resolutions. Further- 
more, whilst errors in determining the intensity from very 
optically thick regions have little effect upon the level popu- 
lations in their neighbourhood, they can strongly affect the 
appearance of an object because they affect the shape of 
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the boundary that can be probed by radiation of a particu- 
lar wavelength. 



3 BENCHMARKING 

In order to verify the accuracy of a code, it is necessary 
to compare its output against either analytical solutions or 
the results of other codes. We tested TORUS against a one- 
dimensional collapsing cloud model to benchmark the iter- 
ative solver used to ascertain the non-LTE level population 
solution along the radius of the cloud. In order to test the 
datacube generation, we chose a two-dimensional, very op- 
tically thick circumstellar disc model in LTE. We compared 



this against results obtained by mcfost (Pinte et al.|2006 



2009) 



3.1 A collapsing cloud in HCO+ 

This problem was first presented as a robust test of a molecu- 
lar line radiative transfer code's ability to deal with as many 
astrophysical phenomena as possible at a 1999 workshop on 
'Radiative Transfer in Molecular Lines' held in Leiden. It 
has since become a standard test of a molecular line trans- 
fer code's ability to reproduce level populations in a com- 
plex physical situation. Although the spherically symmetric 
model is quite straightforward to implement, the velocity 
and temperature gradients, variable turbulent line widths 
and multiple levels provide a stern test of a code's accuracy, 
especially at higher optical depths. It is noted that whilst no 
analytic solution exists for this benchmark, multiple codes 
have reached a broad consensus (within 20 per cent, though 
often much better depending on the fractional abundance of 
HCO"*" and the specific level being verified) as to what the 
level populations should be along the radius of the cloud. 
The results of this test were published in |van Zadelhoff et al.| 
(20021). 



The problem is based on a model by Rawlings et al.| 



around a protostar. The model is similar to the inside-out 
collapse model theorized by |Shu| ( |1977[ ). In this problem, 
the collapse is signified by the negative radial velocity in the 
cloud interior. Because the collapse can only propogate out- 
ward at the sound speed in the cloud, the outer shell of the 
cloud is static (see Figure [I|. The density profile follows a 
piecewise power-law n{r) = no(r/ro)~'" where ro = 10^^ cm 
and m = 1.5 for r < vq and m = 2 beyond tq. All other pa- 
rameters are specified at 50 logarithmically spaced points 
between 10^^ cm and 4.6 x 10^^ cm except the relative abun- 
dance of HCO"*" to H2 which is constant at either [HCO"*"] 
= 10"^ (model (a)) or [HCO+] = 10"** (model (b)). 

TORUS uses an AMR grid which is refined according to 
a series of user-specified rules. In this problem, where the in- 
put parameters are tabulated, this grid is split so that no cell 
contains more than one datum point as defined in the source 
model. In each octal, the salient parameters, temperature, 
n(Il2), systematic velocity field and local turbulent velocity 
are assigned assuming conditions at the centre of the cell 
apply throughout its extent. In this problem both thermal 
and turbulent velocities are combined. Where the cell cen- 
tre does not coincide with a datum point the parameters are 




Radius (cm) 



Radius (cm) 




Radius (cm) 



Radius (cm) 



1 1992 ) to analyse HCO"'' data for an infalling envelope et al 



Figure 1. Parameters as a function of radius along tlie cloud. 
The stark change in n(H2), temperature and velocity at 10^^ cm 
is evidence of the collapse front having reached this radius. 



logarithmically interpolated where an analytical value is not 
available. 

To maintain consistency with other results we used the 
molecular data supplied with the problem set; a database of 
Einstein coefficients and coUisional rate coefficients of HCO"'' 
with H2 in J = ( |Monteiro|[T985l |Green|[T975| . For radu 
inside and outside the region covered by the data file we 
assumed a vacuum at Tcmb. 

The resultant level populations as calculated by TORUS 
are compared for the optically thin scenario (model (a), Fig- 
ure [2) and for the optically thick scenario (model (b). Fig- 
ure Isl) with an average of other codes (see [van Zadelhoff 



2002 for details) that completed this test. The first 



8 levels are all converged to better than 0.25 per cent al- 
lowing us to be confident that the first 6 are accurate. It 
is therefore these first 6 levels that are compared. Higher 
levels are more sparsely populated and are subject to sig- 
nificantly more Monte-Carlo noise. It can be seen that the 
level populations are consistent with other codes. Without 
the exact values obtained for each code it is not possible to 
precisely compare our output with their ensemble average 
but our code varies by no more than 5 per cent for J = in 
the optically thick case. Agreement in the optically thin case 
is even closer. We note that the average of the results of the 
codes may not be a good quantitative measure of accuracy 
but expect that agreement will be closest with other AMC 
codes. 

Using Ng acceleration, the optically thin case converges 
to better than 1 per cent (RMS all levels) in less than 5 
minutes on a single-processor desktop machine (2 GHz Intel 
Core 2 Duo) having reached 12800 random rays using the 
ray-doubling method described in section[2] Within 10 min- 
utes this improves to better than 0.5 per cent using 51200 
rays to sample the radiation field. While each level converges 
at its own rate, in this model the highly populated J — 
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Figure 2. Relative level populations for J = to J = 5 for the optically thin bcnehmark case (tj^2— i ~ 6.5) where [HCO"* 
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Figure 3. Relative level populations for J = to J = 5 for the optically thick benchmark case (tj^2— l ~ 44) where [HCO"^] = 
10~®. Note the deviation from the benchmark of the level populations at the inner boundary reflects differences in the chosen boundary 
conditions of the models that were combined to produce the benchmark curves. 
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Table 1. GG Tau model parameters 



Parameter 


Value 




(unit) 


no 


6.3 X 


109 


cm~^ 


i^CO/Ha 


1.76X 


10-8 




^0 


100 




AU 


^0 


14.55 




AU 


To 


30 




K 


Vo 


3.3 




kms"^ 


''turb 


0.2 




kms~^ 


Inclination 


43 




o 


■'"in 


180 




AU 


rout 


800 




AU 



and J = 1 levels converge much faster than the other levels 
having converged to better than 0.05 per cent after 51200 
rays. 

The relatively low abundance of HCO^ in the model 
creates a scenario where the cloud is quite optically thin. 
The (1-0) transition {vg = 89.1885 GHz) has an optical 
depth of ^ 5 from the centre to the edge of the cloud for 
[HCO"''] — 10~9. The critical density for this transition is 
~ 10^ cm~^. Consequently, throughout most of the cloud 
radiative processes are dominant over coUisional processes 
therefore the assuming LTE would be inappropriate. Fur- 
thermore, the presence of a systematic velocity field allows 
radiation emitted at the line centre of a particular transition 
to escape the denser cloud centre more easily, pushing the 
solution further from LTE. Figure [4] shows the departure co- 
efficients from the LTE case, defined as n"^^^/nf^'^, the 
ratio of the calculated relative level population for a specific 
rotational energy level at some radius to that predicted by 
assuming the cloud is collision dominated. 



3.2 GG Tau - A raytracing test 

The putrey, Guilloteau, fc Simon] |T994t model of GG Tau, 
a young multiple star system surrounded by a circumbinary 
disc, provides an excellent test of the imaging capabilities of 
TORUS. The circumbinary disc is assumed to be in LTE and 
the temperature, H2 number density and velocity profiles 
are all given in analytical form unlike the tabulated data for 
the collapsing cloud in section [3T| This allows a direct com- 
parison of the predictions of the fiux that should be observed 
from the object in ^^CO by different RT codes. We compare 
our results with those of Pinte (private communication). A 
description of the code used to obtain those results is given 



Pinte et al. (20061. The model parameters and scaling 



laws common to both simulations are given below and in 
Table m 



n{r, z) = no 



exp 
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Figure 5. The discretised density profile of GG Tau in g cm~'^. 
The log scale captures the steep decline in material away from 
the disc midplane. 



As the model is symmetric about the rotational axis, a cylin- 
drical coordinate system (r, cf), z) is used. TORUS takes advan- 
tage of rotational symmetry by projecting any coordinate 
with (/!> 7^ onto the (r, 0, z) plane. Naturally, this reduces 
the number of grid cells required to fill a space to 0{h?) 
from 0{h^). In order to resolve the disc well it is necessary 
to devise criteria that ensure the grid is split sufficiently to 
capture all its features. We used the following conditions to 
decide whether to split a cell, where d is the cell width and 
H{r) is the characteristic scaleheight of the disc at radius r: 



H(r 

z 



< 5 and 



d 



H{t) 



>0.1 



^5 and 3 > 5 
H{r) d 

r > 0.99rin and r < l.Olrin and d > O.lnn 

If any of these conditions were true then the cell should 
be split. The cell is not split if r < 0.99rin or r > l.Olnn even 
if another condition is satisfied as the model is not defined 
outside these radii. The minimum cell depth is 3 and the 
maximum depth is set to 20, although this is never reached. 
The first two conditions cover the entire vertical extent of the 
disc and stipulate that at least 50 cells must be used to cover 
the first 5 scaleheights above the disc midplane. Beyond this, 
where the disc is far more tenuous, the criterion is relaxed 
so that the maximum cell size is less than 20 per cent of 
the z— coordinate of the cell. The third criterion ensures the 
adequate resolution of the inner-edge of the disc. 

When the disc is split according to these criteria, a 
cylindrical region with diameter of 6 x lO'^^cm « 2000 AU 
is split into 23,782 cells, the smallest of which is 0.12 AU 
across. Figure [5] shows the discretized density profile and 
Figure l6| shows the disc imaged in ^^CO (1-0). 

Figure W\ illustrates the excellent agreement between 
MCFOST and TORUS line profiles in this test case. The ge- 
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Figure 4. Departure coefficients for 5 different abundances from [HCO^] = 10~^ to 10~^ (see legend). Both model cases [HCO"*"] 
= 10~^ and [HCO+] = 10^* show that they are far from LTE. As expected, for [HCO+] > 10~^ it can be seen that the cloud has 
become sufficiently optically thick in most lines that these levels are well approximated by an LTE solution everywhere except the cloud 
periphery. 




Figure 6. This three colour plot of GG Tau in ^^CO (1-0) shows 
that the observer is inclined to the midplane (as evidenced by 
the elliptical evacuated region in the centre of the disc) and that 
the disc is relatively optically thin almost throughout. The red 
component is material moving away from the obsever, blue to- 
wards and green has no radial velocity component. The intensity 
is given in erg s"'^ cm~^ sr"^ Hz"'^. 



ometry and simplifying assumptions make it possible to test 
the capabilities of the rendering code. Having passed both 
tests we are confident in attempting a far more complex ge- 
ometry that does not exhibit any symmetry and where no 
parameter is defined analytically. 




Velocity {kms"^) 



Figure 7. This line profile shows the excellent agreement be- 
tween MCFOST (solid grey line) and TORUS (dashed black line). 
A distance of 150 parsecs was assumed to convert intensity into 
flux. The agreement is well within observational error and dis- 
crepancies may be put down to minor differences in discretiza- 
tion strategy. This confirms that the TORUS raytracing routines 
are sound. 



4 AN EFFICIENT PARTICLE-MESH 
ALGORITHM 

In order to map the density structure of a particle-based 
representation onto the adaptive mesh employed by TORUS 
it is necessary to use a particle-mesh algorithm that interpo- 
lates data stored at irregularly spaced points onto the grid. 
Many methods exist that already do this; the simplest is to 
use only the particles within a cell to determine a parameter 
value. For example, the density is taken to be the average 
density over all the particles within the cell or the sum of 
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the particle masses divided by the cell volume (Kurosawa 
et al.|[2004 L However, the variation between these two cal- 
culations may be large and neither takes into account the 
concept of smoothing inherent in the SPH technique. Fur- 
thermore, whilst this technique tends to give an accurate 
result in high density regions, where the particle density is 
naturally greatest, it is not capable of giving an answer in 
regions where no particles exist and further provision must 
be made to evaluate parameters in empty cells. More so- 
phisticated algorithms take into account a particle's con- 
tribution to its associated cell and its nearest neighbours 
using a linear or quadratic kernel, namely the cloud-in-cell 
algorithm, however even this algorithm fails to take into ac- 
count the variable radius over which each individual particle 
contributes to its environs. In fact, in order to accurately de- 
termine information about a parameter at any given point 
in space from a particle ensemble, it is necessary, in theory, 
to calculate and sum the contribution from every particle. 
In practice, to make SPH calculations tractable, each par- 
ticle is assigned a characteristic size scale that defines its 
'sphere of influence', beyond which the contribution is pre- 
sumed to be negligible. This variable smoothing length, hi, 
for a particle with index i, is either determined by 

1/" 



hi ^ r] 



(13) 



where 77 is a constant controlling the approximate number 
of nearest neighbours a particle has and v is the dimension- 
ality of the simulation ( [Springel &c Hernquist|2002 Price fc 
Monaghan 2004, 2 007[ ), or is determined at each timestep to 
control exactly the number of nearest neighbours ( Bate et al 



2003 1, depending on the desired resolution of the calculation. 



Note that in the latter method, no closed form exists for hi. 
The weight of each particle's contribution is determined by 
the smoothing kernel; a normalized analytical function of 
the particle-point separation in units of smoothing lengths, 
WdP — pi\, hi). Throughout this paper we use the cubic 



spline kernel (see equation 151. This kernel is popular due 
to its ease of computation and compact support over 2h. 
Kernels are discussed in greater detail in Monaghan] ( |T992[ ). 
Below we outline the method we use to map the information 
stored in the particle representation onto the AMR grid. 

4.1 Creating the grid 

An AMR grid is created by the repeated subdivision of an 
initial unrefined cell centred on the entire region of interest. 
The parent is bisected once in each dimension, thus each 
child has a volume 2"" times that of its parent. As we are 
dealing with complex astrophysical structures in this section 
we use a full 3-dimensional representation of the space so 
each parent has 8 child cells, or octals. More detail of the 



AMR method is given in Harries (i2000 ) and Symington et al, 
(20051). 



Problems in star formation typically span many orders 
of magnitude in both space and density. SPH and AMR 
naturally resolve the fine structure of these problems well 
whereas a fixed, regular grid will often be insufficient to 
capture crucial information like peaks in the density profile 
where a protostellar core far smaller than the size of a grid 
cell has begun to accrete material. The flexibility of AMR is 
that it will adapt the resolution given to a region depending 



on the criteria the user applies. Typically, one is interested 
in resolving the temperature and density profile in a region, 
although when calculating line radiative transfer accurately 
it is necessary to accurately resolve variations in the velocity 
field as well owing to the anisotropy in the absorption intro- 
duced by the line profile function. To create the grid we test 
each parent cell against a mass criterion, a density criterion 
and a velocity criterion. Where these criteria are not met, a 
cell is recursively split until the conditions in each child cell 
are met up to a specified maximum depth of recursion. 

The mass per cell criterion determines the maximum 
number of particles that may occupy a cell, thus ensuring 
that no cell has more than a certain fraction of the entire 
mass of the simulation within it. Typically, this value is set 
to ~ 50 particles; this balances the need for accurate deter- 
mination of densities within cells and sufficient grid resolu- 
tion within memory constraints. 

The density criterion dictates the range of particle den- 
sities within a cell. This criterion facilitates the resolution 
of density gradients which are both physically and compu- 
tationally important in radiative transfer calculations (e.g. 
treatment of shocks, jets, etc.). Large changes in density 
from cell to cell can have deleterious effects on the con- 
vergence of the iterative algorithm used to determine the 
relative level populations of a molecular species within the 
cell so it is crucial that the volume over which they change 
is well resolved. Typically, this criterion is set so that the 
maximum density stored on a particle within the cell is not 
greater than twice that of the minimum density. 

The velocity criterion regulates the variation in the 
magnitude of the velocity within a cell. Again, it is vital to 
have good velocity resolution over a cell to obtain the cor- 
rect line profile shape since absorption and emission of line 
radiation is a strong function of Doppler-shifted frequency. 
Therefore, where possible we use this criterion to ensure that 
the range of velocities within a cell is less than 5 times that 
of the turbulent line velocity (see equation [7|. 

If any of the conditions in the parent cell meet the split- 
ting criteria then the cell is split. This process is iterated 
until each cell no longer triggers any of the criteria or a 
maximum cell depth is reached. This is to ensure that the 
simulation does not run out of memory. We impose a further 
condition that the minimum subcell size is no smaller than 
the smallest inter-particle separation (~ 0.25^min), which is 
commensurate with the natural resolution of the SPH sim- 
ulation. Figure [S] shows a cut-through version of a mesh 
created using the above splitting criteria. The minimum cell 
size in this grid is ~ 0.05 AU. 

4.2 Determining cell parameters 

Having created a grid it is then necessary to populate each 
cell with parameters pertinent to the calculation that is to 
be performed. The method outlined here is optimized to take 
advantage of geometry where possible but is not prejudiced 
against the general case where the distribution of particles 
is not known a prion. 

For each cell in the grid, the centre is chosen to be the 
point in space where the conditions that pervade the cell 
are to be determined. Whilst it is trivial to determine the 
region over which a particle acts (because the distance at 
which a particle ceases to contribute is fixed at 2hi), the 
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Figure 8. A cut-through of an adaptively refined grid. The min- 
imum AMR depth of 3 can clearly be seen here giving a coarse 
grid of 8^ cells; the deepest level of recursion is 23. The spatial 
extent of this grid is 5 X lO"*^* cm. 



inverse problem of determining which particles contribute 
to the grid cell centre is less trivial. That is, some may lie 
within the cell but some may not. Furthermore, some may 
lie within the cell yet not contribute to it at all. Naively, one 
might assume that by finding the particle with the greatest 
smoothing length, /imax and subsequently checking for all 
particles that lie within 2/iinax of a point one is guaranteed 
to recover the most accurate answer possible. Whilst this is 
certainly true, it turns out to be so inefficient that deter- 
mining the contents of an entire grid like this would take 
an impractically long time. Not only does /imax necessarily 
overestimate the smoothing lengths of all but one particle 
in the ensemble, but also, due to the direct relationship be- 
tween smoothing length and density it is unrepresentative 
of the entire physical region; very low density volumes be- 
ing the exception rather than the rule in circumstellar discs 
and star-forming clusters. Moreover, these very low density 
regions play a negligible role in radiative transfer and are 
naturally undersampled by the SPH technique. We define 
the parameter /i95% as the smoothing length of the particle 
lying at the 95"^ percentile of the distribution of smoothing 
lengths of all particles. This length is often as much as an or- 
der of magnitude less than that of /imax but still greater than 
that of all but the 5 per cent of particles representing the 
least well-populated areas of the physical space. However, 
as the adaptive mesh is already split so that the cell size is 
commensurate with the smoothing lengths of the particles it 
contains, a better estimate for the critical lengthscale over 
which to search for nearby particles is rcrit = min(4d, ^195%), 
where d is the cell width. The factor of 4 has been empir- 
ically derived as the smallest factor that retains the fur- 
thestmost particles that are likely to contribute whilst keep- 
ing computational speed high. Equation [13] shows that for 
each ten- fold reduction in r, 3 orders of magnitude of den- 



sity resolution are lost. However, as the radius over which 
contributing particles are searched for reduces, the volume 
(and hence the computational effort) reduces as r^, making 
the problem more accessible. Furthermore, as it is the least 
dense regions that are no longer resolved, the impact on the 
solution is minimised. If however there are no particles in 
a point's vicinity then the code makes another attempt to 
locate a nearby particle that can be used to derive the lo- 
cal density. The algorithm is repeated with a wider search 
volume whose radius is determined by 



min(max(4d, 2/195%), 0.5/in 



(14) 



This second condition covers almost all cases where the ini- 
tial condition does not suffice but if even this radius is not 
sufficient then one last attempt is made using rcrit = /imax; 
beyond this the cell is justifiably declared empty and its 
parameters set to some global minimum. 

Sorting the particle list by ascending a;— coordinate fa- 
cilitates the reduction of the number of candidates that 
can contribute to a point, P. This is achieved by removing 
all particles that lie outside 2r-crit in the abscissal dimen- 
sion, leaving only particles whose x— coordinate lies within 
Pi ± ^crit. In the case where all the particles are equally 
spaced this reduces the number of comparisons necessary to 
4 ^/Ai'part . Each particle in this reduced list is then tested 
against the same criterion in y, i.e. |P^ — pf | < 2rcrit. For 
all those particles where this criterion is satisfied a further 
test for the z— component is made. This process discovers all 
particles that lie within a cube of side length 4rcrit centred 
about the point. Whilst this means that only approximately 
half the particles are expected to lie within a sphere of ra- 
dius 2rcrit , it is still more efficient than directly calculating 
the pairwise distance for each particle in the entire ensem- 
ble. The logic is such that, at each step in the algorithm, 
it is necessary to do less work than the stage before. Fur- 
thermore, the order in which the criteria are applied has 
been selected in such a way as to reject the greatest number 
of particles as early on in the process as possible; that is, 
culling by x then y then z is the most efficient method of 
locating nearby particles for a disc with its semi-major axis 
in the xy-p\a.ne. Once the reduced list of particles has been 
created it is then a matter of calculating the contribution of 
each particle to the point. 

For each remaining particle, the contribution to the 
point is determined by the kernel smoothing function, W . A 
common smoothing kernel is the cubic spline, defined below. 



W{q,h)^ 



where 



W(g) 
/i" 






1 - |g^ 4 




for ^ q < 1, 

for 1 ^ q < 2, 

otherwise, 

(15) 



Qi = 



(16) 



|P-PJ 
hi 

and o is the normalization constant for the dimensionless 
kernel W(g) defined such that j^ W{q)dV = 1. For 1/ = 3, 
a = I/tt. By utilizing the property of compact support pos- 
sessed by the kernel smoothing function it is possible to dis- 
card those particles with q ^ 2 without having to evaluate 
the kernel function. 

In a homogeneous, isotropic distribution of 10® parti- 
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cles, this smoothing kernel reduces the number of particles 
able to contribute to a point by approximately five orders of 
magnitude. Specifically, the number of particles able to con- 
tribute is independent of volume or number, being constant 
at ~ 60 for r/ = 1.2. 

In the SPH method, the value of a scalar parameter at 
a point in space, A(r), is the weighted sum of the stored 
parameter values, Ai, over all nearby particles, i.e. 



^(r) = E 



Pi/i" 



W{q,)Ai 



J?-" J2 w(gO^« 



(17) 



by equation 1 13| (18) 



This formulation assumes that X]W(qi) ~ 1 which is typ- 
ically true for r in the bulk of the simulation surrounded 
by a full complement of neighbours. Where it is true it is 
possible to normalize A{r) by dividing by ^ W((?i) ~ 1, re- 
ducing the variation associated with kernel smoothing. The 
superposition of many spherically bound kernels on to a reg- 
ular 3D grid will not recover a constant field but rather it 
will tend to oscillate around the true value; the field will be 
over-estimated at the particle position (where the weighting 
function at that point is 1 and the sum over all other parti- 
cles is positive) and will similarly be underestimated at the 
point in space centred at the midpoint of surrounding ker- 
nels. This normalization can certainly improve the appear- 
ance of an image, (see iPrice, 2007) but can cause problems 
at free surfaces where a smooth decay is preferred. 

Particles with h > ^/Kim/A^part, the average linear sep- 
aration of particles in the simulation volume, Kiim, are ex- 
pected to lie close to a free surface and are initially classed 
as 'hull' particles. The algorithm adaptively decides whether 
to normalize or not depending upon one of two strategies: 
either, if '^W{qi) > 0.3 then normalize; or normalize only 
if all contributing particles are determined to be 'bulk' par- 
ticles, (as opposed to hull particles). Either strategy relies 
upon empirically determined parameters, '^W{qi) > 0.3 
typically recovers total masses most similar to the SPH sim- 
ulation. The designation of a hull particle may change if it 
is subsequently found to belong to a set of particles where 
X] W(gi) > 0.3, whereby it is classed again as a bulk particle. 
Similarly, it is possible, though unlikely, that a bulk particle 
may become a hull particle if it belongs to a set that fails 
to meet the criterion. In any case, it is possible for a user to 
override this strategy as it is not obviously desirable to let 
the temperature or velocity decay at a free surface even in 
the apparent absence of matter. 

Once the contributing particles have been located and 
their weights determined they can be reused for each pa- 
rameter that needs to be found. A vector such as velocity 
must be split into its 3 components and each one must be 
determined before outputting the resultant vector. The al- 
gorithm is repeated for each point until the grid has been 
filled with data. 



5 IMAGING SIMULATED STAR FORMATION 
REGIONS 

We have taken snapshots o f an SPH star formation simu- 
lation by [Bate et"aL] (|2003[). The simulation used 3.5 x 10** 



particles to trace the evolution of a 50 Mq molecular cloud. 
The cloud was initially spherically symmetric with a uni- 
form density and temperature (10 K). A supersonic turbu- 
lent velocity field {M — 6.4) was imposed on the cloud to 
create anisotropies which ultimately led to the formation of 
a dense pre-stellar core after approx imately 1 free-fall time. 
Here, ts ^ 1.9 x 10^ 



yr- 



Bate et al 



observed that high den- 
sity regions form from converging gas flows during the decay 
of the turbulent velocity field. These enhancements continue 
to increase in density until a particle exceeds a density of 
10~^^g cm~^, when a sink particle is created with a mass to 
the sum of the masses of all particles within a 5 AU radius. 
The sink particle inherits a weighted average of all the prop- 
erties of the contributing particles and they are destroyed. 
This ameliorates the computational slowdown caused by re- 
solving the acceleration caused by large forces generated by 
particles with small separations which requires a very small 
timestep. 

During the early stages of collapse, the relatively tenu- 
ous gas allows energy to be radiated with ease. At this stage 
the gas was assumed to be isothermal. If the rate at which 
energy is released exceeds the rate at which the gas can cool 
then the local environment begins to heat rapidly. [Masunaga] 



& Inutsuka ( 2000 1 showed that the gas begins to heat sig- 
nificantly at a density of p > 10"^"^ g cm~^ depending on 
gas opacity (and initial temperature). To avoid performing 
a computationally expensive full radiative treatment, |Bate| 
|et al.| used a barotropic equation of state p — Kp^ . The 
barotropic exponent, rj, equals unity where the gas is isother- 
mal (for p ^ 10~^^ g cm~^) increases to 7/5 for p > 10"^^ 
g cm~^. This equates to a temperature-density relationship 
of 



T 



max 10, 10 



vio-i«/ 



2/5 



(19) 



The frames we have used in this paper are taken at 1.0, 
1.1, 1.2, 1.3 and 1.4 ta; the first snapshot occurs just before 
the start of star formation, the second is taken during the 
first star formation episode, while the third and forth are 
taken during an interstitial period during which a second 
and third star-forming core is formed. The final snapshot is 
taken at the end of the simulation when over 50 objects have 
been formed, some of which are still accreting. At the end of 
the simulation, approximately 1/5 of the SPH gas particles 
(« 10 Mq) have been accreted onto the sink particles. 



5.1 Transposing SPH particles onto the AMR 
grid 

We discretized each SPH snapshot using the algorithm out- 
lined in section |4] using the same splitting criteria for each 
timestep; no more than 50 particles were allowed within one 
cell, the ratio of the maximum density to the minimum den- 
sity must be no greater than 2 and the range of velocity 
magnitudes must be less than 5uturb- Figure [9] illustrates 
how an SPH representation of density can be discretized 
onto a mesh. 

The process of discretization will necessarily introduce 
some further uncertainty in the exact distribution of mass 
throughout the cloud, however the algorithm has been de- 
signed to minimize any systematic error which may arise 
through using too coarse a grid to represent the space, 
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Table 2. Table of SPH properties and corresponding AMR grid properties for each timestep. 
SPH properties Grid properties 



Timestep 



No 



N„ 



pmax 



Noctals Max Depth Mass (Mq) Mass Error (per cent) Npart/Noctals 
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Figure 10. Histogram of the relative fraction of AMR cells (dark 
grey) at each depth for i = l.ltff. The smoothing lengths of the 
SPH particles (light grey) were binned according to their closest 
association in log-space. Note that the x-axis is reversed to show 
the proportion of cells at lower depths on the left. 



Figure 9. A slice through the xy— plane illustrating the dynamic 
range of densities that can be captured within the cloud using 
AMR. Thirty orders of magnitude are present in the simulation 
with densities reaching as high as 2 X IQ- ^g cm-'^ and as low 
as 5 X IQ- ^^g cm-'^ (although very little material by mass is 



Our analysis of this molecular cloud seeks to exam- 



actually present below > 
plot is given in n(H2). 



10" 



'). The greyscale for this 



a deeper study of which is undertaken in |Acrenian et al.| 
1 2010b I. The algorithm has been successfully used on galax- 



ies (Acreman et al. 2010a I, circumstellar discs ( Acreman 
|et al.||2010b[ ), star-forming clusters (this work) and as a 
crude test of discretization accuracy the total mass stored 
on the grid is found to lie within 5 per cent of the true mass 
in all cases. In this work, the mass discretization error is 
better than 3 per cent and in the absence of overdense re- 
gions, better than 1 per cent (see Table [2|. Furthermore, 
a comparison of the number of cells in each level of refine- 
ment against a histogram of the smoothing lengths of the 



particles (Figure 10 1 shows broad agreement, except in the 
highly refined tail of the AMR cell distribution. This may be 
because of the additional constraints required for adequate 
radiative transfer, additional velocity splitting in particular. 
Over all 5 frames, the maximum level of refinement was 22, 
resulting in the smallest octal having a linear dimension of 
^ 0.08 AU, commensurate with the shortest inter-particle 
spacing. On average, this discretization policy resulted in a 
ratio of particlesxells of approximately 1/3. 



ine the conclusions of Ayliffe et al. (20071 following the 



analysis of Walsh et al. (20041. Accordingly we use the 
i^CO (1-0) transition (z/q = 110.2013542798 GHz) to trace 
low density envelopes. Initially, we assume an abundance of 
[I'^CO] — 1Q~^ relative to II2. To trace gas closer in to the 
cores within the cluster we use the C'^^O (1-0) transition 
{uo = 109.7821734 GHz, [C^^^O] = lO-'^) and to define cores 
we use the N2H+ (1-0) transition {vo = 93.1737 GHz) be- 
cause of its high critical density (ncrit ~ 1.4x 10^). For N2H"'' 
we assume a constant relative abundance of 1.5 x lO-^"^ 
in line with values obtained for well-studied cores in the 
Taurus- Auriga molecular complex ( [Tafalla et al.|[2004[ ). For 
all molecules, we assume a non-thermal turbulence parame- 
ter with a fuUwidth of 0.3 km s-^ ~ Cs for H2 at the simu- 
lated temperatures. 

We superimposed the ^^CO, C^^O and N2H+ molecular 
data on to the discretized grid for a total of 15 distinct 
grids (3 molecules and 5 timesteps). We derived a non-LTE 
solution for each grid using the molecular line transfer code 
outlined in section [2] Table |3] shows the mean error of the 
least well converged level populations (J = 4 in all cases) 
over all grid cells. A further condition that every cell had 
to have an RMS error of less than 1 per cent in J = and 
J = 1 was imposed to reduce any pixel-to-pixel variance in 
intensity in the subsequent raytracing of the grid. 

To improve the statistics, each of the 15 grids were ob- 
served from 20 equally spaced positions on a sphere with 
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Table 3. Fractional error (in per cent) of the least well converged 
level in each grid for each timestep. 

Species 





13CO 


CiSQ 


N2H+ 


1.0 tff 


0.30 


0.14 


1.03 


1.1 tff 


0.27 


0.19 


1.98 


1.2 tff 


0.25 


0.07 


0.54 


1.3 tff 


0.16 


0.12 


0.59 


1.4 tff 


0.14 


0.15 


0.71 



radius 140 pc. Care was taken to maximize the separation of 
each observer position giving equal coverage of the viewing 
sphere from each perspective. A 1024 x 1024x80 element dat- 
acube was created for each observer position-tracer-timestep 
triplet. This work therefore contains 300 distinct datasets 
upon which to perform statistical analysis. The maps have 
a linear spatial resolution of 2 x 10^^ cm or ~ 133AU per ele- 
ment. This is equivalent to a spatial resolution of ~ 1 arcsec 
at 140pc. From the known velocity distribution of the SPH 
particles, we expect emission to come from a range of veloc- 
ities within 3.2 km s~^ of the lin e c entre giving a velocity 
resolution of 0.08 km s" 



Figure 11 shows an example im- 
age that can be produced by combining the non-LTE line 
transfer algorithm, the SPH particle interpolation and the 
raytracing routine. In these grids, molecular abundance has 
been assumed to remain constant in both time and space. 
Without a detailed evolution history, it is not possible to pre- 
dict with any certainty the local chemistry in any region of 
the grid. In section [5^ we describe a simple chemical model 
to take into account the expected CO freezeout in the cloud. 
In the absence of better chemical modelling, we assume that 
the abundance of N2H^ remains constant, which is appro- 
priate given the age of the cluster ([Bergin &: Langer|1997|). 



5.2 Chemistry 

Molecular adsorption onto dust grain surfaces is known to 
be an effective depletion mechanism in many car bonaceous 
species at densities of n(H2) = 3 x 10* cm~^ | Bacmann 
et al.|[2002 K Notable absentees in this list are N-containing 
species including N2II"'". C^*0 is expected to trace a larger 
scale than N2H+ not only because of its lower critical density 
but also because it effectively freezes out at densities above 
~ 3 x 10* cm~^. Therefore, in the cluster, C^*0 intensity is 
weighted towards the outer regions while N2H"'" traces the 
denser inner parts. 



5.2.1 Drop profile 

To test the effect of CO chemistry we use a 'drop' pro- 

(20041). In its sim- 



J0rgensen et al. 



file as described in 
plest implementation, the model assumes a constant, un- 
depleted abundance of CO relative to that of H2, Xo, where 
Tdust/gas ^ Tovap Or n(H2) < 3 X 10* cm~^ combined with a 
depleted abundance of Xo, elsewhere. Using equation ( |19[ ), 
it is possible to reformulate the model solely in terms of the 




Figure 11. A cluster observed in HCO+ (1-0)- Blue represents 
material moving towards from the observer; green, stationary ma- 
terial and red receding gas. The smooth transitions in intensity 
are due to the careful choice of normalization strategy, adaptive 
pixel sampling and density subsampling. 



local gas density, p, 



X{p) 



Xo 
Xu 



for 10"^^ < p< 1.56 X 10 
elsewhere, 



'^ gem ^ 



(20) 

Empirically, the drop profile has been used to improve the 
fit of Gaussian line profiles to low J observational C^*0 
data of various class 0, class 1 and pre-stellar objects. The 
improvement in the goodness of fit is found to be sufficiently 
significant as to provide strong evidence for its acceptance 
over the constant abundance model. Moreover, the model 



has its basis in observations ( Caselli et al. 1999 Tafalla et al 



2002 1 where a freeze-out zone has been directly imaged. 



Physically, the profile is explained in terms of radii. The 
material within the inner radius is expected to be gas-phase 
owing to its proximity to the warm core. The canonical tem- 
perature at which CO is thought to have returned to the gas- 
phase is 30 K. Beyond the outer radius, at densities lower 
than 3 x 10* cm~^, the depletion timescale due to freezeout 
is thought to exceed the lifetime of the core and is thus as- 
sumed to remain at an undepleted abundance, to reduce the 
number of free parameters. 

Since the critical densities of N2H+ are higher than 
those of the CO isotopologues, they are less sensitive to the 
outer region of the envelope where depletion and contribu- 
tion from the surrounding cloud may be important. It has 
been posited ( J0rgensen et al.]|2004 1 that this may explain 
why these transitions can be modelled assuming a constant 
density. 

In our chemistry runs we repeat the same method as 
outlined in section |5.1| but with depleted densities one or- 
der of magnitude lower in the intermediate region, viz. 
[1^CO]d = 1 X 10-'' and [CI'^OId = 1 x 10"^ 
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5.3 Analysis 

For each synthetic observation the intensity datacubes had 
their background intensity subtracted, _B(i-'o,rcmb), and 
were summed over all velocity channels to give an integrated 
intensity map. The maps were converted to brightness tem- 
perature, Tb using the Rayleigh- Jeans approximation, i.e. 



Ti 



2ku2 



(21) 



We used wavdetect, part of the CIAO data analysis 
system written for the Chandra X-ray Observatory, to de- 
tect cores in the N2H+ intensity maps. The wavdetect tool 
( [Freeman et al.|2002| | correlates the image with wavelets of 
different, user-defined scales and then searches the results for 
significant correlations. This method is very similar to the 
difference of Gaussians (DoG) method, where a difference 
map is created by convolving an image with two Gaussian 
distributions with different spatial variances and subtract- 
ing the resultants, preserving structure with lengthscales be- 
tween the two variances, wavdetect produced very similar 
results but was slightly better at detecting objects that it 
was not possible to resolve using DoG and was more robust 
against the converse case of rejecting objects that upon later 
inspection had distinct line profiles (often indicative of con- 
taining protostars). We looked for cores with lengthscales of 
2-64 pixels (~ 0.001 — 0.04 pc) and applied a threshold 
of 50 per cent of the global maximum to the transformed 
image to determine the membership of a pixel to a core. A 
recursive blob counting algorithm was invoked to determine 
the number of separate regions identified as containing a lo- 
cal maximum within the map; for N2H^, each region is a 
potential core and for '^^CO and C^^O each region defines 
the envelope which surrounds the core. This is similar to 
the method of |Offner, Klein, fc McKee| ( |2008a[ ) in that we 
generate a background map and then detect local intensity 
peaks above that. 

Any core candidate that lay wholly or partially outside 
the enveloping material as determined by the low density 
contour was excluded from the analysis. In this constant 
abundance analysis, it is expected that the regions should 
be strict subsets of each other. For each core candidate, a 
mask was applied to the original datacube corresponding 
to the region defined by the 50 per cent contour and line 
profiles for that region were found by integrating over all 
pixels in the region for each molecular tracer. 

The line profile was fitted with a Gaussian distribution 
using a Levenberg-Marquardt algorithm. The initial Gaus- 
sian was set to have a peak brightness temperature of 5 K, 
o" = 1 km s"'^ and to be centred at u = km s~^. Each da- 
tum point was assumed to have equal weighting. In addition 
to characterizing the line profile using those parameters, the 
mean and median of the distribution were calculated pro- 
viding a robust method to verify the position of the peak in 
the Gaussian distribution determined by fitting. Similarly, 
the full-width at half-maximum (FWHM) of the profile was 
recorded to parameterize the variance. The same process 
was repeated for the low density tracer and the resultant 
difference in line centre velocities was assumed to represent 
the systematic shift in line-of-sight velocity. Figures |12[ |13| 
and [14] illustrate the results of the analysis for one particular 
observation at 1.4 tfi. 



5.4 Results 

Having removed the two cores that did not lie within the 
50 per cent low-density envelope we also removed 11 cores 
(over the 5 timesteps) which were judged to be spurious de- 
tections which were typically characterized by very small ar- 
eas (< 10 pixels), similar velocity profiles to a larger nearby 
detection and often did not contain any sink particles. This 
left 169 cores over the 5 timesteps; the fewest detected cores 
(22 from 20 observing positions) being in the first timestep 
before the creation of any sink particles and the most (58) 
being detected in the last timestep, where three unique star- 
forming regions were detected by Bate et al. (20031. Every 



remaining valid core-envelope pair was analyzed and the re- 
sults were collated to allow us to perform a statistical anal- 
ysis according to the strategy outlined in section [53] Cores 
that contained one or more sink particles were judged to be 
protostellar cores while cores without sink particles are cat- 



egorized as being prestellar (or starless). Following Ayliffe 
et al. (20071, we do not present the results of Gaussian fits 



in this work. All analysis was done using both methods and 
the differences between methods of analysis were sufficiently 
small that they did not affect any conclusions drawn upon 
the results of the work. Moreover, the profiles of ^^CO data 
were rarely found to be Gaussian, exhibiting flat, wide peaks 
(owing to their large optical depth). As such, we character- 
ize the dispersion using the FWHM / 2\/21n(2) of any given 
line profile and take the line centre to be the median, m, of 
the distribution (J™ f{x) dx = 0.5 J^ f{x) dx). As a non- 
parametric estimator of the line centre, the median does not 
assume normally distributed data; as a result the median is 
less influenced by shoulders, long tails, noise etc. and may 
give a less biased estimate of the true line centre. Conse- 
quently we did not attempt to flt multiple Gaussians to the 
data (with the assumption that a double-peaked distribu- 
tion indicated two cores along the line-of-sight) primarily 
because it was possible to confirm from the SPH particle 
distribution that this was not the case but also because it 
allows us to more directly compare our results with pre- 
vious studies. The presence of a doubly (or more) peaked 
profiles is caused by the velocity shift associated with the 
rotation and/or infall of gas. In order to illustrate these ef- 
fects we have calculated profiles from 25 areal bins around 
the cores shown in Figure [12] and plotted them over an inte- 
grated intensity map (Figure [T5| . Whilst the optically thin 
gas is associated with single-peaked profiles, the denser cores 
display a more complex line profile morphology. A similar 
variety of profile shapes have been found in one-dimensional 



simulations of collapsing and rotating cores ( Pavlyuchenkov 
et al. 12008] ). 



Figure [Tm shows a histogram of the distribution of rela- 
tive velocity difference between N2H+ cores and low density 
gas envelopes identified by ^^CO summed over all timesteps. 
The standard deviation of all cores (both pre- and proto- 
stellar) is a = 0.16 km s~^. In the same figure, we plot 
three normalized Gaussian distributions whose standard de- 
viations are equal to that of the equivalent standard devia- 
tions as derived from the mean FWHM of the line profiles 
for each tracer molecule. They are notably more broad than 
either the relative velocity difference or the sound speed. 

Figure [17] illustrates values for the mean velocity dis- 
persions of each tracer for each timestep. The trend for the 
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Figure 12. Integrated intensity map in N2H+ (1-0). Contours have been added to demarcate each of the cores identified. The fine 
profiles of the two largest cores and those of the enveloping gas are illustrated in Figures |13| and |14| respectively. Note that this image 
only shows the region of interest in the cloud. The extent of the entire cloud is ^ 0.6 pc. 
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Figure 13. Typical core profiles. Dashed lines indicate the Gaussian profile fit and FWHM. Dotted line indicates position of the Gaussian 
line centre. The legend also displays median line centre and the equivalent standard deviation, Ce, from the FWHM of the profile, [(b) | 
shows a profile that it would be possible to fit two Gaussians to (assuming two cores along the line-of-sight) however by examining the 
SPH particle distribution this profile is known to be caused by rotating gas surrounding a multiple system. 



CO isotopologues is the same; that the mean velocity dis- 
persion is relatively low at 1.0 ts, rises until a peak at 1.2 
tff and then turns over and levels off towards the end of 
the simulation. The same effect can be seen in the N2H"'" 
cores except that the effect seems to be delayed until 1.3 tff. 



We discuss this trend in section [S] The relative differences 
in line-centre velocities between the C^*0 envelope and the 
N2H+ cores are lower throughout the cluster. This is due 
to C^*0 emission being more optically thin than the ^^CO 
emission. Consequently, C^*0 (1-0) is able to trace a region 
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Figure 14. Typical low-density envelope profiles traced using ^^CO. Both profiles exhibit flat peaks associated with the saturation of 
optically thick ^^CO (1-0) line. 
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Figure 15. The line profiles of the square regions are superimposed onto a N2H+ (1-0) integrated intensity map (see Figure IT2I. The 
regions containing the three cores have somewhat broader line profiles and exhibit more structure (i.e. double peaks caused by rotation 
and infall) than the optically thin gas surrounding them. 
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Figure 16. Histogram of velocity difference of all cores (combined 
starless and protostellar cores) identified at all time steps (light 
grey, a = 0.16 km s~^) and protostellar cores (dark grey, a = 0.18 
km s~^). Also plotted are the velocity dispersions for N2H+ (dot, 
o- = 0.861 km s-l), Cl*0 (dot-dash, a = 0.808 km s"!), ^^CO 
(dash, a = 1.05 km s"-"^). 
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Figure 17. The velocity dispersion of each molecular species for 
each timestep. White squares denote N2H+, grey squares denote 
^■^CO and black squares denote C^^O. 



closer to the core. Note that in these constant abundance 
calculations CO is not frozen-out and therefore will trace a 
higher density region. 



5.4-1 The effect of chemistry 

We repeated the analysis with chemical model outlined in 
section |5.2| Although minor differences in line profile in- 
tensity and shape are observed, the velocity dispersions did 
not change significantly and no systematic differences with 
the constant abundance models were detected. It appears 
that the optical depths of the transitions we have studied 
are such that the effective photo-surface of the cores is out- 
side the volume where freeze-out occurs. Naturally a more 
sophisticated treatment of chemistry, involving a full time- 
dependent calculation, is required to confirm our prelimi- 
nary results. However such a calculation is well beyond the 
scope of this paper. 



6 DISCUSSION 

We performed a detailed radiative transfer analysis of a self- 
gravitating hydrodynamics simulation in order to examine 
the relationship between the core and envelope velocity dis- 
persions and make a more direct comparison between the 
hydrodynamical models and millimetre observations. We 
used information about specific molecular tracers and de- 
termined the non-LTE level populations and used these to 
calculate the emergent intensity from different observation 
angles around the cluster. We found good qualitative agree- 
ment with observational results (e.g. [Walsh et aL|2004[|2007| 
KirketalJ2007) , our main conclusion being that one cannot 



reject competitive accretion as a viable theory of star forma- 
tion based on observed velocity profiles. ) Ay liffe et al.| ( [2007| 
found that the majority of their sources had velocity differ- 
ences (between high and low density gas) less than the sound 
speed at all timesteps but that they had a significant tail out 
to large velocity differences (in this case > 0.5 km s~^). How- 
ever, they also found that their high-density gas linewidth 
became larger than their low-density gas linewidth at times 
greater than 1.1 tfj, contradicting the observations of|Walshl 



et al.land Kirk, Johnstone, & Tafalla (2007*). We believe this 



may bo because the study done by Ayliffe et al. (20071 did 



not take into account optical depth effects and that their 
density cut for the high density material was too high. By 
reducing their assumed critical density by a factor of 3 they 
were able to reduce the velocity dispersion of high-density 
material by a significant amount so that the high-density gas 
linewidth was smaller than the low-density linewidth, when 
averaged over all timesteps. This work does not show the 
trend of larger high-density linewidths at any time, although 
Figure [17] shows that the velocity dispersion turns over with 
increasing time. We believe that this effect can be explained 
by two opposing processes. The original simulation did not 
drive turbulence so it is expected that in the absence of dy- 
namic interactions between protostellar objects the velocity 
dispersion of the gas will ultimately tend to decrease. How- 
ever, as we are observing regions undergoing complex inter- 
actions that act to stir up the gas, thereby increasing the 
velocity dispersion, we see the two effects combined; each 
star formation episode acting to increase the velocity dis- 
persion and decay over time acting to reduce it. Like | Ayliffe| 
|et al.| and the observational surveys, we see no clear trend 
in the relative motions of core/envelope pairs and note that 
the standard deviations of the velocity differences are always 
far smaller than the velocity dispersion. 

The original hydrodynamic simulation did not include 
any form of driven turbulence. We have assumed a global, 
constant, non-turbulent line width of 0.3 km s~^. This is 
similar to that of the mean non-thermal turbulent line width 
of the driven case in [Offner et al.| ( |2008b[ ) but at odds with 
the theory that cores are supported by thermal pressure be- 
cause large non-thermal motions will have disappeared on 
small scales (< 0.1 pc) and at high densities (> 10^ cm" 



(Johnstone et al.|2010l. As a result we found larger velocity 



dispersions than have been observed in some studies ( Offner 
|et al.||2008aj . However, not only are our observed core mo- 
tions smaller than the mean velocity dispersion but also 



they are often less than the local sound speed (Figure 16 1. 
This is true for both starless and protostellar cores and in 
fact, our results mirror those of|Kirk et al.|in that they too 
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show no indication that the core-to-envelope motions signifi- 
cantly change between the starless and protostellar stages of 
evolution. By performing a more detailed radiative transfer 
study of |Bate et al.[ s cluster simulation, our findings go some 
way towards reconciling the results previously obtained by 



Ayliffe et al. (20071 with those found by others who had 



performed similar studies. 

Although relatively sophisticated we intend to incorpo- 
rate additional physics in future. Specifically w e have not yet 
taken into account N2H''" hyperfine splitting. Tafalla et al. 



(20041 studied the effect of neglecting hyperfine splitting 
(which maximizes radiative trapping). They found that the 
difference in obtained level population was on the order of 
tens of per cent which is comparable to the uncertainties 
in collision parameters. Neglecting hyperfine splitting can 
have a significant effect on the optical depth of any N2H"'" 
transition and consequently the emergent intensity. In ef- 
fect, we have systematically overestimated the line optical 
depth and the inclusion of hyperfine splitting would shift 
the line formation regions to higher density. In the future 
we wish to incorporate hyperfine splitting, although previ- 
ous work ( [Daniel et al.|2006i ) incorporating this microphysics 
has demonstrated that whilst it is likely to have a significant 
effect it will not be large enough to affect our overall con- 
clusions. 

Another important factor is depletion of CO by ad- 
sorption onto dust grain surfaces. Throughout this paper, 
we have assumed that all tracer molecules occur with con- 
stant abundance and many observational studies have shown 
that this is demonstrably not the case. Chemical networks 
(Bergin fc Langcr 1997', [Glover fc Mac Low 2007) exist that 



allow us to predict the abundance of many molecular species 
in molecular clouds. We used a simple chemical model to 
estimate the effect that CO chemistry would have on our 
observations, and have discovered that a simple freeze-out 
model does not significantly affect our results. Nonetheless 
we recognize that a more complete treatment is desirable 
and a natural extension of this work would be to couple 
TORUS with a chemical network to understand more fully 
where line emission comes from in molecular clouds. 



throughout the cluster for tracers of high- and low-density 
gas over 5 timesteps covering the creation of 50 protostellar 
objects we were able to produce and analyse synthetic obser- 
vations. We have shown in this paper that clusters exhibiting 
competitive accretion are able to reproduce the properties 
of relative core motions found in observation. 

As well as enabling us to accurately model line pro- 
files, the mapping of density distributions from the particle- 
based representation onto an AMR grid has important im- 
plications for coupling radiation transfer and hydrodynam- 
ics. For example, using the technique outlined in section |4[ 
it is now possible to conduct a full, polychromatic multiple 
scattering treatment of the transport of radiation in an SPH 
simulated circumstellar disc using high-accuracy grid-based 
techniques ( Acreman et al.|2010b \ without having to resort 
to the flux-limited diffusion approximation employed in ear- 



lier radiative transfer simulations (e.g. Whitehouse & Bate 



2004[ [Krumholz et al.|2007[[Bate[[2009b[ ). However, extend 
ing such calculations to a cluster collapse model is currently 
not tractable given the very large dynamic range of linear 
scales necessary to accurately treat the radiation transport 
of optically thin-to-thick boundaries, but, with further ad- 
vances in computational power this kind of calculation will 
become a reality, bringing with it enhanced certainty in the 
conclusions that can be reached using these simulations. 
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7 CONCLUSION 

This paper demonstrates the molecular line radiative trans- 
fer capabilities of TORUS, a radiative transfer code that uses 
AMR to resolve the fine details of complex three-dimensional 
environments. The molecular line transfer module was com- 
pared against other codes in a scenario where LTE condi- 
tions did not apply and the results showed excellent agree- 
ment. 

In environments such as molecular clouds, gas that is 
dense enough to emit with sufficient intensity in a given 
molecular line may only exist in small regions when com- 
pared with the size of the whole cloud. Consequently, fine 
discretization is only required over a small volume. Thus, we 
devised an efficient method of mapping irregular SPH data 
onto an AMR grid with minimal discretization error. This 
technique permitted the study of the expected molecular 
line emission from a hydrodynamical simulation of molecu- 
lar cloud collapse using SPH with our grid-based radiative 
transfer code. Having obtained non-LTE level populations 
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